/* slfprmpf.cpp by K.Tsuru */
// function ID 4108 DRADIX since ver 2.182
// first version made on 30 Apr 2010
/**************************************************************
 permutation
          n!
 nPr = ---------
       (n - r)!
is evaluated by factorizing n! and (n - r)! into prime factors
**************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
// for sorting
static int sortByPower(const void *e1, const void *e2) {
    return (int)((const primeFactor *)e2)->power - (int)((const primeFactor *)e1)->power;
}

SLong permL(const ulong n, const ulong r) {
    if(r == 0) return 1.0;
    if(n == r) return Fact(n);

    if( n < r ) SNManager::SetError(SNManager::DOMAIN_ERR,"combPF()", 4108);

    SNBlock<primeFactor> npf, n_rpf;
    int in   = FactorialIntoPrimeFactor(n, npf);        // number of prime factor n!
    int in_r = FactorialIntoPrimeFactor(n - r, n_rpf);// (n-r)!

    // reduction
    for(int k = 0; k < in_r; k++) npf[k].power -= n_rpf(k).power;
    int non0pwr = 0;    // number of non zero power elements
    for(int k = 0; k < in; k++) {
        if(npf(k).power) non0pwr++;
    }
    // pick out nonzero power elements
    primeFactor* wpf = new primeFactor[non0pwr+1];
    int el = 0;
    for(int k = 0; k < in; k++) {
        if(npf(k).power) {  wpf[el] = npf(k);   el++;   }
    }
    wpf[el].prime = wpf[el].power = 0;  // index of last element

    qsort(wpf, non0pwr, sizeof(primeFactor), sortByPower);

    SNBlock<primeFactor> pfBk(el+1);

    for(int k = 0; k <= el; k++) pfBk[k] = wpf[k];

    SLong p = PrimeFactorProduct(pfBk, el);
    delete[] wpf;
    return p;
}